use rand::Rng;

fn gcd(a: u128, b: u128) -> u128 {
    if b == 0 {
        a
    } else {
        gcd(b, a % b)
    }
}

fn modmul(a: u128, b: u128, m: u128) -> u128 {
    let mut res = 0;
    let mut a = a % m;
    let mut b = b;
    while b != 0 {
        if b & 1 == 1 {
            res = (res + a) % m;
        }
        a = (a * 2) % m;
        b >>= 1;
    }
    res
}

fn modpow(x: u128, p: u128, mod_: u128) -> u128 {
    let mut ans = 1;
    let mut x = x;
    let mut p = p;
    while p != 0 {
        if p & 1 == 1 {
            ans = modmul(ans, x, mod_);
        }
        x = modmul(x, x, mod_);
        p >>= 1;
    }
    ans
}

// Miller-Rabin primality test
// Returns true if p is probably prime, false if it's definitely composite
fn miller_rabin(p: u128) -> bool {
    if p < 2 {
        return false;
    }
    if p == 2 || p == 3 {
        return true;
    }
    let mut d = p - 1;
    let mut r = 0;
    while d & 1 == 0 {
        r += 1;
        d >>= 1;
    }
    for _ in 0..10 {
        let a = rand::thread_rng().gen_range(2..p - 1);
        let mut x = modpow(a, d, p);
        if x == 1 || x == p - 1 {
            continue;
        }
        for _ in 0..r - 1 {
            x = modmul(x, x, p);
            if x == p - 1 {
                break;
            }
        }
        if x != p - 1 {
            return false;
        }
    }
    true
}

fn pollard_rho(x: u128) -> u128 {
    let mut s = 0;
    let mut t = 0;
    let c = rand::thread_rng().gen_range(1..x - 1);
    let mut step = 0;
    let mut goal = 1;
    let mut val = 1;
    loop {
        for _ in 0..goal {
            t = modmul(t, t, x) + c;
            if t >= x {
                t -= x;
            }
            val = modmul(val, (t as i128 - s as i128).abs() as u128, x);
            if step % 127 == 0 {
                let d = gcd(val, x);
                if d > 1 {
                    return d;
                }
            }
            step += 1;
        }
        let d = gcd(val, x);
        if d > 1 {
            return d;
        }
        goal *= 2;
        s = t;
        val = 1;
    }
}

fn factor(x: u128, max_factor: u128) -> u128 {
    let mut x = x;
    if x <= max_factor || x < 2 {
        return max_factor;
    }
    if miller_rabin(x) {
        return max_factor.max(x);
    }
    let mut p = x;
    while p >= x {
        p = pollard_rho(x);
    }
    while x % p == 0 {
        x /= p;
    }
    let max_factor = max_factor.max(factor(x, max_factor));
    max_factor.max(factor(p, max_factor))
}

pub fn find_max_prime_factor(number: u128) -> u128 {
    factor(number, 0)
}
